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ABSTRACT 

In the standard model of cosmic structure formation, dark matter haloes form by 
gravitational instability. The process is hierarchical: smaller systems collapse earlier, 
and later merge to form larger haloes. The galaxy clusters, hosted by the largest dark 
matter haloes, are at the top of this hierarchy representing the largest as well as the last 
structures formed in the universe, while the smaller and first haloes are those Earth- 
sized dark subhaloes which have been both predicted by theoretical considerations and 
found in numerical simulations, though it does not exist any observational hints of their 
existence. The probability that a halo of mass m at redshift z will be part of a larger 
halo of mass M at the present time can be described in the frame of the extended Press 
& Schecter theory making use of the progenitor (conditional) mass function. Using the 
progenitor mass function we calculate analytically, at redshift zero, the distribution of 
subhaloes in mass, formation epoch and rarity of the peak of the density field at the 
formation epoch. That is done for a Milky Way-size system, assuming both a spherical 
and an ellipsoidal collapse model. Our calculation assumes that small progenitors do 
not lose mass due to dynamical processes after entering the parent halo, and that they 
do not interact with other subhaloes. For a ACDM power spectrum we obtain a subhalo 
mass function dn/dm proportional to m~ Q with a model-independent a ~ 2. Assuming 
the dark matter is a weakly interacting massive particle, the inferred distributions is 
used to test the feasibility of an indirect detection in the 7-rays energy band of such 
a population of subhaloes with a GLAST-like satellite. 
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1 INTRODUCTION 

The present-day description of the universe includes the 
presence of a large amount of cold dark matter (CDM) whose 
nature and distribution is unknown. This Dark Matter (DM) 
provides about 26 % of the energy budget of the universe. 

The amount and properties of CDM is well constrained 
by astrophysical observations such as the anisotropies in 
the Cosmic Microwave Background, l arge scale structure 
and distant typ e I A supernovae ()Spergel et al.l 120031 : 
lAstier et al1l2006l : iTegmark et alll2006j ). On the other hand, 
two main open questions arise. The first concerns the par- 
ticle physics nature of the CDM. Weakly interacting mas- 
sive particles (WIMPs) are attractive ca ndidates since thei r 
relic abundance can fit the observed one [|Dimopouloslll990h . 
Stable neutralinos in supersymmetric extensions of the stan- 
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dard model (SUSY) jjungman et all Il99rj ; | Bertone et alj 
l2005al ) or Kaluza-Klein particles (KKP) in theories with 
a TeV" 1 size universal ex tra dimension (|Appelquist et all 
l200ll : IServant fc TaitJliool) are the most commonly studied 
particles. Since these particles have never been observed, 
there is a large uncertainty on the prediction of their effects 
which has to be taken into account. The other open question 
regards the distributio n of DM inside the haloes. Numeri- 
cal N-body simulation s (N avarro et al.lll997l :l Diemand et al.l 
l2004d : iNavarro et al.ll200ot ). whose scale resolution is about 
~ 0.1 kpc, allow solely an extrapolation of the very inner 
slope of the DM profile and do not take into account in- 
teractions with the baryons which f all in the DM potentia l 
well or the presence of inner cores (jBerezinskv et al.l 120031 ) 
or the controverse effect of the presence of a black hole 
at the centre of the halo dUllio et al. 200 ll: Bertone et alj 



l2005bl : lBertone fc Merrittjl2005l : lMerritt et alJl2002T ). Expert 
imental data on DM distribution in the haloes of galaxies 
and clusters are not conclusive too (see, i.g., the discus- 
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sion in iFornengo etU] l|2004 )). In the hierarchical forma- 
tion scheme of the CDM scenario, large systems are the 
result of the merging and accretion of smaller haloes (sub- 
haloes), whose dense central cores would survive the merg- 
ing event and continue to orbit within the pa rent halo, as 
shown by high resolution N-body simulations (jMoore et al.l 
ll999l ; lGhigna et al.ll200ol ; iBlasi fc Shethll2000l ). CDM mod- 
els are characterized by an excess of power on small scales. 
The arising divergence of the linear density constrast at 
large wavenumbers has been proved to be damped by colli- 
sional processes and free streaming, respectively before and 
after kinetic decoupling, leading to exponential damping 
of the linear CDM density contrast and to the existence 
of a typical scale (Jeans scale) for the fir st haloes corre- 
sponding to a Jeans mass ab out 10 _6 Mq (Ho fmann et al.l 
l200ll ; iGreen et all |2004 120051 ). Numerical simulations have 
indeed reproduced hierarchical clustering in CDM cosmolo- 
gies with a mass resolution sufficient to resolve the Jeans 
mass l|Diemand et alj|2005al ) with particle mass m p = 1.2 x 
10 -10 Mq and force resolution of e = 0.01 pc; however such 
a high resolution run could be evolved only to z = 26, in 
a very small spatial patch, and producing haloes of mass 
[1O- 6 ,1O- 4 ]M . 

Among the simulations evolved on larger scales and to 
redshif t z = 0, present mile stones are the Millennium Simu- 
lation (|Springel et al.ll2005l ) and the Via Lactea Simulation 
l|Diemand et al. 2007a! ) . The first is a cosmological N-Body 



run with over 10 billion particles in a cubic region 500 Mpc/h 
on a side (particle mass m p = 1.23x 1O 9 M0; force resolution 
e = 7 kpc) ; the second was done to obtain a simulated Milky 
Way with the highest possible mass resolution (particle mass 
m p = 2.09 x 1O 3 M0; force resolution e = 90 pc). However 
a simul ation with the mass a nd force resolution similar to 
that of l|Diemand et al . 2005a), evolved to redshift zero over 
a region containing a mass comparable to that of our Galaxy 
would require about 10 20 particles and a time resolution of a 
few years. Such requirements are way beyond the computa- 
tional capabilities of present-day supercomputers: applying 
Moore's law and starting from present day state-of-the-art, 
a run like this could be performed in roughly 50 years from 
now. 

A reasonable alternative is to study the clustering prop- 
erties of Milky Way-like systems through an analytical ap- 
proach. We use the fact that the probability that a halo of 
mass m at redshift z will be part of a larger halo of mass 
M at the present time is described by the progenitor con- 
ditional mass function f(m,z\M, zo — 0), according to the 
so-called extended Press & Schechter theory. Using the pro- 
genitor mass function, we can calculate analytically, at red- 
shift zero, the distribution of subhaloes in mass, formation 
epoch and rarity of the peak of the density field at the for- 
mation epoch. That is done for a Milky Way-size system, 
assuming both a spherical and an ellipsoidal collapse model. 

Nu merical simulations described in iDiemand et alJ 
(2005b) show that the distribution of material originating 
from the earliest branches of the merger tree within the 
present day haloes depends on the a-peaks of the primordial 
density fluctuation field it belonged to. We extend their nu- 
merical results by performing an analytical estimate of the 
density peaks distribution as a function of the halo mass 
traced back to the smallest scale haloes, thus avoiding the 
limitation imposed by numerical simulations. In this way 



we obtain a realistic estimate of the distribution and mass 
function of the whole population of subhaloes. 

Such an analytical estimate can provide a powerful tool 
to take into account the effect of early high-density peaks in 
present day haloes. 

This is particularly important in the framework of dark 
matter indirect detection, since a high cr-peak halo translates 
into a higher concentration and thus a higher value for the 
density squared which has to be integrated along the line of 
sight to obtain a prediction for particle fluxes coming from 
dark matter annihilation. 

Given some model for the hierarchical formation of our 
Galaxy, and for the internal structure of subhaloes, DM may 
be in fact indirectly detecte d using annihilation rates pre- 
dicted from particle physics l|Bergstromll200ol ; iBertone et al.l 
l2005al ) through the observation of high density point-source 
or extended regions inside our Galaxy. If we restrict our- 
selves to 7-ray observations, these can be obtained using ei- 
ther atmospheric Cerenkov telescope s ([Weekes et al 1 ll997l : 
lAharonian et al.lll997|: | Ba ixcras 2003) or satellite-borne de- 
tectors like GLAST (|Morselli et al.lll997l ). The detectability 
of DM substructures with GLAST has been w idely discussed 
in the literature (see, e.g. IPieri et all l)2007h and references 
therein). The small mass haloes have been found to give the 
main contribution to an unresolved 7-ray foreground aris- 
ing from DM annihilation, while their detection as resolved 
objects has been proved to be very unlike. Indeed the un- 
resolved subhalo foreground is prominent above the MW 
smooth foreground far from the Galactic Center, where the 
overall flux is still too low to be detected. 

In this paper we apply the analytical derivation of the 
subhalo population properties, such as the cr-peak distribu- 
tion, on the indirect detection of 7-rays. We thus study the 
possibility that high cr-peak material could arise the fore- 
ground level above the detectability threshold of a GLAST- 
like large field of view satel lite. 

As in IPieri et al.1 (|2007l ). we use different models for the 
virial concentration of subhaloes. 

The paper is organized as follows: in Sec. O we review 
the spherical and ellipsoidal collapse model and their prop- 
erties. In Sec.|3]we describe the original analytical derivation 
of the density peak distribution as function of the halo mass, 
and the subhaloes mass function for a present day halo with 
mass M = 10 12 Mq/H. In Sec. [4] we estimate the upper 
bound for the contribution to the 7-ray flux due to the pres- 
ence of a population of subhaloes inside the Milky Way. In 
Sec. we study the prospects for detection of substructures 
with a GLAST-like experiments in our best case scenario. A 
discussion of our results can be found in Sec. [S] 



2 EXTENDED-PRESS & SCHECHTER: FROM 
PROGENITORS TO SUBHALOES 

In the hierarchical picture of galaxy formation, structures 
up to protogalactic scale grow as a consequence of re- 
peated merging events. Smaller systems collapse at high red- 
shifts, when the universe is denser, an d subsequently assem - 
ble to form bigger and bigger haloes |Lacev fc Cole! ll993T ) . 
This merging history is often represented by the so called 
" merger-trees" . 

Smaller systems accreted onto a larger halo along its 



© 2007 RAS, MNRAS OOO.rTlfTol 



Analytical Subhaloes Population 3 



merging-history-tree and still surviving at a later time are 
called "substructures" or "subhaloes" ([Ghigna et al. 199S ; 



Tormen et ai1l2004l; 



van den Bosch et al 



Gao et allliooilDe Lucia et al.M2004a 



20051 ). In what follows we will discuss 



an analytical approach to derive the mass function of sub- 
haloes. We will use the simplifying assumption that no tidal 
stripping nor merging events among substructures happen. 
In this approach the mass of each subhalo r emains constant 
in time, and equals the original virial mass (|Eke et al.lll996l ) 
of the progenitor halo at the co n sidere d redshift. A similar 
study was carried out by IShethl ([2003), who calculated the 
subhalo mass function using the creation rate of the pro- 
genitors of a present day dark matter halo; our approach 
is different: we derive the subhalo mass function from the 
entire population of progenitors (as shown by Eq. [4} , in or- 
der to allow a d i rect co mparison with the JV-Body results of 
iDiemand et"ail l|2005bl ). 



2.1 Conditional Mass Function 

Let us consider a halo with virial mass M at some final 
redshift zo- According to the hierarchical picture of galaxy 
formation, going backward in time the halo will be splitted 
in smaller and smaller systems, called "progenitors". Mass 
conservation tells us that the sum of all masses of progenitor 
haloes at any given redshift equals the mass of the halo at Zq. 
Let us define the conditional mass function f(m, z\M, zo)dm 
as the fraction of mass belonging to haloes with mass be- 
tween m and m + dm at redshift z, which are progenitors of 
a halo of mass M (a M-halo) at a later redshift zq. 

Assuming the spherical collapse model 
|Press fe Schechterl Il974f ), we can express m and z as 
a function of the new variables s and 5 SC . The conditional 
mass function is independent on the powe r spectrum of 
densi ty fluctuations and it is described as (|Lacev fc Cold 
1 19931 ): 



f(s,S 3C \S, S )ds 



exp • 



(8 SC — 6p) 2 
2(s - S) 



ds 



where s ■ 
a m-halo 



W 

= a 2 (rn) is the square of the mass variance of 
and S ac is the spherical collapse overdensity 
at redshift z. S and (5o are the mass variance of an M- 
halo and the spherical collapse overdensity at the present 
time, respectively. To compute the mass variance we have 
chosen a power spectrum with primordial spectral index 
n = 1, and a transfer func tion obtained from CMBFAST 
ijSeliak fc Zaldarriagalfl996h for a concordance ACDM uni- 
verse (f2 m , f2 A , h = 0.3, 0.7, 0.7) with <t 8 = 0.772, extended 
down to a mass M = lO 6 M /h. 

We have integrated this power spectrum using a top- 
hat filter in real space. To obtain the mass variance until 
the typical Jeans neutralino mass we linearly extrapolate 
the log(m)-s relation to M = 10~ 6 M e /h. 

Over the last ten years TV-Body simulations have shown 
that the collapse of dark matter haloes is actually not well 
described by an isolated spherical model; the influence of 
sourrounding proto-haloes can be reproduced using an ellip- 
soida l model jSheth. Mo fc TormenfeoOll : ISheth fc Tormenl 
120021 ). 

In the excursion set approach, the progenitor mass func- 



spherical collapse 
-ellipsoidal collapse 



z = 0.5 




-10 
log(m/M) 

Figure 1. Ellipsoidal (solid) and spherical (dotted) conditional 
mass function computed for a present-day dark matter halo with 
mass 10 12 Mq//i and for five different redshifts. 



tion of a halo is described by the conditional probability of 
first upcrossing distribution. Such a probability is well fitted 
by a random walk in the plane (s, 5), starting from (S, 80) 
l|Bond et al.lll99ll ). 

In the spherical collapse model this barrier has a con- 
stant height, defined by the collapse redshift: B 3C (s,S sc ) = 
5 SC . For the ellipsoidal collapse case the barrier height is not 
constant, but depends on s and on S sc as described by the 
following equation: 

B ec (s,5 sc ) = ^«5 sc [l + / 3(-^-) 7 ] . (2) 

ISheth. Mo fc Tormenl (|200ll ) found q = 0.707, /3 = 0.5 and 
7 = 0.6; the value of the last two parameters is motivated 
by an analysis of the collapse of homogeneous ellipsoids, 
whereas the value of q comes from requiring that the pre- 
dicted halo abundances match what is found in the simula- 
tions. 

Considering the barrier described in Eq[2] 
ISheth fc Tormenl (2002) found an approximate solution for 
the diffusion equation, expressed as follows: 



f(s,5 ac \S,S )ds = 



\T(a,6 ac \S,6 )\ x 
>/27r(a - S) 

[B{s,S sc ) 



(3) 



exp • 



ds 



2(s - S) 



with T(s\S): 
T(s,8 sc \S, So) 



= E 



(s- 



s) n d n [B(s,5 s 



B(S,S )] 



ds n 



In Figure [T] we show the conditional mass function at 
five different redshifts for a halo with present-day mass 
M = 10 12 Mq/H, both for the spherical (dotted curves) and 
ellipsoidal (solid) collapse prediction. It can be observed that 
the halo is splitted in smaller and smaller progenitors at 
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higher redshifts; discrepancies between the two models de- 
pend both on mass and on redshift. 

Comparing the two prediction at fixed redshift, one 
can note that the spherical model predicts more progen- 
itors at intermediate mass, and fewer at both very small 
and very large masses, compared to the ellipsoidal model 
l|Sheth fc Tormenll2002h . In other words, the two predictions 
cross each others in two points, although these crossings do 
not necessarily fall in the range of masses plotted in the 
Figure. 

A direct consequence of this is that massive progenitors 
exist at higher redshifts in the ellipsoidal collapse, and the 
distribution of formation redshifts (defined as the earliest 
epoch when a halo assembles half of its final m ass in one sys- 
tem) is consequently shifted to earlier epochs (|Giocoli et al.l 
l2007h . 

From f(s, S SC \S, 5o)ds we can write the total number of 
progenitors at any given redshift as: 

N{m,5 sc \M,S )dm = ^lf( s ,5 sc \S,S a )ds. (4) 
m(s) 

Considering a scale free power spectrum P(k) oc k n , the 
mass variance scales as s(m) oc m - '™" 1 " 3 ^ 3 , and the number 
of progenitors can be explicitely written in terms of s: 




N(m,S BC \M,S )dm = 



(n+3)/3 



f(s,S BC \S,S )ds. (5) 



Figure 2. Total number of progenitors in a given mass bin, 
as a function of redshifts, for a present day halo with mass 
M = 10 12 M®/h. For each mass bin we show the prediction for 
spherical (dotted lines) and ellipsoidal (solid) collapse models. 



2.2 Number of progenitors 

Interating Eq. [3] over mass we obtain the total number of 
progenitors in the given mass interval, as a function of red- 
shifts: 



dn(z, Am) 



r m f 

J rm 



N{m,S sc \M,S )dm = N(z) 



(6) 



where rrii and mj represent the bounds of the interval. For 
a white-noise power spectrum (scale free with n — 0) and a 
spherical collapse mass function, a primitive of this integral 
can be written as: 

, {Ssc - S ) 2 

N(z) = — l —\e 2(s-5) ( 7) 



(S sc ~S ) 2 



SV2^y 



[2V^~S(Ssc - So) - e 2(s-S) 



2n[S - (S sc - So) lerf 



8 



V^(s - S) 



)] 



In Figure [5] we show the total number of progenitors 
in five different mass decades, for a halo with mass M = 
10 12 M©//iat 20, as a function of redshifts. We have assumed 
a concordance ACDM power spectrum and have integrated 
Eq.[4]numerically. The solid lines represent the prediction for 
the ellipsoidal collapse model while the solid ones refer to the 
ellipsoidal model. From top to bottom the curves represent 
the following mass bins: [h 1CT 6 , 1(T 5 ], [10 _1 ,1], [10 2 ,10 3 ], 
[10 6 , 10 7 ] and [10 9 , 10 10 ], all but the first expressed in term 
of M @ /h. 

It can be observed that the spherical collapse, for a 
fixed mass bin, underpredicts the number of haloes at high 
redshifts compared to the ellipsoidal model. We will see in 



the next sections that if we consider the variable v{z, m) = 
S S c(z)/a(m), for any given mass this will result in the in- 
equality u ec (m) > v sc (m). 



3 UNEVOLVED SUBHALOES MASS 

FUNCTION FROM THE MERGER TREE OF 
A PARENT A/-HALO 

The progenitors mass function, integrated over S sc , gives the 
total number of progenitors of mass between m and m + dm 
that a halo of final mass M has had at all times: 



dn(m) 
dm 



/ — f{s,8 sc \S,5o)dS s 

JSn m 



(8) 



in the case of the spherical collapse this integral results in: 



dn(m) 
dln(m) 



M [da/dm] 
2tt yf7~S~ 



(9) 



with a as I for a LCDM power spectrum. Since the same sys- 
tem may be a progenitor of the same final halo at more than 
one redshift, integrating the progenitor mass function over- 
counts the total number of progenitors. The result of this 
integration must then be prope rly re-normalized by im pos- 
ing the constrain coming from (iDiemand et al.ll2005al ) that 
roughly 10% of the total Milky Way mass (M = 10 12 M@/h) 
is in systems with mass ranging from 10 7 to 10 10 M@/fc 



m A 
M 



0.1 



(10) 



In Figure [3] we plot the differential mass distribution 
of subhaloes in a 10 12 M®/h (Milky Way-like) dark matter 
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Figure 3. Differential distribution of subhaloes in a 10 12 Mq//i 
dark matter halo. The distribution has a slope approximatively 
equal to 1 and has been normalized considering that 10% of the 
total mass is in subhaloes with mass from 10 7 to 10 10 Mq/Zi. 

halo. The distribution has a power law behaviour approxi- 
mately described by the relation: 

^pl = Am-\ (11) 
dm 

with 7 as 2 for both the spherical and the ellipsoidal col- 
lapse model, respectiveljQ. Once fixed the normalization fac- 
tor, we find that the differential distribution of the sub- 
haloes is independent on the mass of the progenitor halo, 
M, considering all the progenitors with mass from 10 -6 Mq 
to m/M = 0.01. 

3.1 Progenitors cr-peak in the host halo 

Using h igh resolution TV-Body simulations, iDiemand et al.l 
|2005bh studied the spatial distribution - at z = - of mat- 
ter belonging to high redshift progenitors of a given system. 
They found that this distribution mainly depends on the 
rareness of the density peak corresponding to the progeni- 
tor, expressed in terms of v = 8 sc /cr(M, z), and is largely 
independent on the particular value of z and M: matter 
from high v progenitors ends up at smaller distances from 
the center of the final system. 

We can understand this in term of the revise d secondary 
infall |Quinn fc Zurekll 19881 ; fZaroubi et aTlll996l '): the forma- 
tion of haloes in TV-Body simulations preserves ranking of 
particle binding energy, that is, particles in the cores of pro- 
genitor haloes will end up in the core of the final system. 
Equally, particles from progenitors accreted at earlier times, 
hence possessing more negative initial binding energies, will 

A least-squares fit on the points gives 7 SC = —1.9972 ± 0.0001 
and 7 CC = -1.9937 ± 0.0003. 



t — i — i — rz 




Figure 4. Progenitor mass function integrated over all redshifts. 
In the top panel we show the distribution for all the masses, while 
in the bottom panel we consider only progenitors in the first and 
last subhalo mass decades. 



likely have a more negative final energy, and so be more 
centrally concentrated than average matter. 

At fixed redshift (hence at fixed S ac ), higher mass pro- 
genitors have a larger v, are more self-bound than smaller 
mass ones, and thus end up closer to the center of the fi- 
nal system. Analogously, for a fixed progenitor mass, higher 
redshift progenitors have a larger 8 SC , hence a larger v\ since 
at higher redshift the universe is denser, they also are more 
self-bound than lower redshift siblings, and so end up closer 
to the center of final system. 

In Figure |4] we plot the subhalo mass function in terms 
of v. To compute the factor v for each progenitor we inte- 
grated the total number of progenitors in a given mass bin 
(Eq. [fjl, at all redshifts. In the top panel we consider all 
the progenitors at all redshifts, with mass in the full range 
/ilO -6 to 10 10 Mq//j; in the bottom panel we show the sim- 
ilar distribution only for the smallest and larger decade of 
the progenitors mass. 
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4 7- RAY FLUX FROM GALACTIC 
SUBSTRUCTURES 

4.1 Modeling Galactic halo and substructures 

We model the distrib ution of DM in our Galaxy after 
iDiemand etafl (|2005bh . 

For the smooth component of the Milky Way we use the 
best fit to the high r esolution numerical experiments of 
IDiemand etaH |2005al ): 



Px( r ) = 



(12) 



with (7, /3,a) = (1.2,3,1). The scale radius r a and density 
p s , are constrained by the viri al properties of the halo. Fol- 
lowing |Dieman^et^al] |2005bT ) we adopt r a = 26kpc, while 
p s has to be normalized to the virial mass of the smooth 
DM halo. We include a physical cutff r cut = 10~ 8 kpc which 
represents the distance at which the self-annihilation rate 
equals the dynamical time of spike formation. 

We shape the spatial distribution of subhaloes accord- 
ing to the fact that it traces the mass distribution of the 
parent halo from r v i r down to a minimum radius r m j n (M) 
where tidal effects become important. We use Eq. [12] to- 
gether with the fact that the dependence from the initial 
conditions when the haloes accreted onto the present-day 
Milky Way halo is set through the p arameter u(M). We the n 
use the parametrization obtained in IDiemand et"al] (l2005bl ): 

T 's * T*v — Ji/Ts 



f v = exp(i>/2) 
— y V = 3 + 0.26V 1 - 6 



(13) 



This parametrization reflects the fact that material accreted 
in areas with high density fluctuations is more concentrated 
toward the centre of the galaxy, and has a steeper outer 
slope. We also use the mass function derived in Sec. [2] to 
model the number density of subhaloes per unit mass at a 
distance r from the GC, for a given v(M): 



p 3 h{M,r,v) 



AM- 2 e(r-r mm (M)) 

a-i (jH u —7) /a ' 



r v {M) 



AM) 



(14) 



Few constraints exist on t he density profile of each sub- 
halo. Numerical simulations (jDiemand et alj l2005al . 120061 . 
l2007bl ) suggest they were formed with a NFW profile, which 
is described by Eq. [T2]with (7, 0, a) = (1, 3, 2). Even if sub- 
haloes probably underwent tidal stripping and consequent 
mass loss after merging, their higher central de nsity should 
preven t the inner regions from being affected. IPieri et"al] 
explored different possibilities for the concentration 
parameter c v i r — r V i r /r 3 , where r v i r is defined as the ra- 
dius at which the mean halo density is 200 times the crit- 
ical density. Following their guidelines, we use two models 
for the concentration c V i r : we assume that the inner struc- 
ture of subhaloes is either fixed at the time they merge onto 
the parent halo (^-labeled model) or that it evolves with 
redshift until the present time (0 model). In model B re /,o 
the NFW concentra tion is computed at z — according to 
iBullock et alj (|200ll ) (hence the prefix B), and extrapolated 
to low masses. In model B re / jZ , the values of c v i r (M,z) are 
obtained from those at z = using the evolutionary rela- 
tion c v ir(M, z) = c V i r (M, z = 0) /(l + z) , where the merging 
redshift z is determined by the knowledge of the value of v 
assigned to each progenitor. Therefore, subhaloes are much 
denser in model B re / jZ than in model B^^o- 

The values c V i r thus found refer to progenitors formed 
from average density fluctuations (u — la peaks of the fluc- 
tuation density field). However, haloes with equal mass at 
redshift z\ may have assembled at different previous epochs; 
specifically, if we call Zi > Z\ the redshift of mass assem- 
bly for progenitors observed at redshift z%, the amplitude 
of the initial density fluctuations producing the progenitors 
is an increasing function of Zi. Therefore, their concentra- 
tion c(M, z) is also an increasing function of the peak am- 
plitude v. To account for this effect, we use the relation 
c v ir{M,v) = u(M)c V j r (M,u = 1), which ha s been tested 
against simulations by IDiemand et all (|2005bh . 



4.2 Modeling the 7-ray flux from Dark Matter 
annihilation 

We model the photon flux from neutralino ann ihilation in 
the p opulation of galactic subhhaloes following IPieri et"all 
(2007). Given a direction of observation defined by the 
angle-of-view ip from the Galactic Center, and a detector 
with angular resolution 9, the 7-ray flux can be parametrized 



= ^(Sy) x $ cosmo (<M) (15) 



dE~, 



The particle physics dependence in Eq. \TE\ is given by the 
annihilation spectrum and DM properties and is embedded 
in the term: 



in units of M^kpc -3 . The mass dependence in r v de- 
pends reflects the mass dependence of the virial parameter 
T s = r V i r /c V i r . The effect of tidal disruption is taken into 
account through the step function 8(r — r m in(M)), where 
r m in(M) is estimated following the Roche criterion. A is a 
normalization factor obtained by imposing that 10% of the 
MW mass is distrib uted in subhaloes with masses in the 
range 10 7 - 1O 1O M l|Diemand et al.ll2005a[ ) as in Sec|2 

As a result about 50% of the Milky Way mass 
is contained within ~ 2 x 10 16 subhaloes in the mass 
range [10 -6 , 1O 1O ]M0. The solar neighborhood density is 
~ 280 pc~ 3 , mainly constituted by haloes with mass of m x is the DM particle mass, a ann u is the self-annihilation 
1O _6 M0. The halo closest to the Earth is expected to be cross-section times the relative velocity of the two annihilat- 
ing particles, and dN^/dE-y is the differential photon spec- 
trum for a gi ven final state / with b ranching ratio Bf, which 
we take from iFornengo etahl (|2004h . 

The line-of-sight integral defined as: 



dE- 



-(£7) 



1 ffann" 

4-7T 2m| 



^dEl Bf - 



(16) 



located ~ 9.5 x 10~ pc away. 

The remaining 50% of the Milky Way mass is assumed 
to be smoothly distributed, and we use this half mass value 
to normalize p s in Eq. [T5] 
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\ip,AQ) = f dM ( dv f [ d6d<f) f dX f dc 

JM Jv J J An Jl.o.s Jc 



[p sh (M,R(R @ ,X,tp,9,<j>),v) x P{v{M)) x P(c(M))x 

x r(A, A', V, 0', A i/, c) x J(x, y, z\X, 6, 4>)] (17) 

accounts for the influence of cosmology in the flux computa- 
tion. AQ. is the solid angle defined by the angular resolution 
of the instrument, J(x, y, z\ X, AQ) is the Jacobian deter- 
minant, R — \/A 2 + Rq 2 — 2XRqC, is the galactocentric 
distance and r is the radial distance inside the single sub- 
halo. Rq is the distance of the Sun from the galactic cen- 
ter and C = cos(6>) cos(V') - cos(^) sin(0) sin(^). P(v(M)) 
is the probability distribution function for the peak rar- 
ity v(M) calculated using the extended Press- Schechter for- 
malism. P(c(M)) is the lognormal probability distribution 
for c centered on c V i r (M) as it is computed in our models. 
While P(u(M) is determined by the merging history of each 
subhalo, P(c(M)) de scribes the scatter in concentration fo r 
haloes of equal mass (|Bullock et al.ll200ll : lNeto et al.ll2007h : 
therefore the two probabilities may be assumed independent. 
The single halo contribution to the total flux is given by 



alo (M, r, v, c) = 



// 

J J A 



d<j>'d0' 

Afi Jl.O.S 



dX' 



Px(M, r(X, X', ip, 0' (/>'), u, 
A 2 



-J{x,y,z\\',0'<t>) 



(18) 



This equation is also used to derive the contribution of the 
smooth component of the MW itself. 

Eq. [T7] gives the average subhalo contribution to the 
Galactic annihilation flux within AQ along the direction ip. 

This contribution is shown in Fig. [5] together with the 
MW smooth halo component obtained with Eq |18l for the 
two models considered in this analysis, for AQ = 10~ J sr, 
corresponding to an experimental angular resolution of 0.1°. 
The sum of the MW smooth and clumpy diffuse contribu- 
tions is shown as well. We define this sum as our "annihila- 
tion signal", which will be multiplied by Eq |16l to obtain the 
predicted 7-ray diffuse flux from neutralino annihilation in 
our Galaxy. In the small box we show a zoom at small angles 
of the an nihilation signal a nd we superimpose the signal ob- 
tained in lPieri et alj ((20071 ) for two similar models (we refer 
to their paper for the detailed explanation of models). Our 
models give a higher flux at the Galactic Center, where the 
signal is dominated by the MW smooth contribution. This 
is due to the different MW profile adopted. Yet, we find 
one order of magnitude of enhancement at the GC in the 
subdominant subhalo contribution as well, due to the pres- 
ence of P(y(M)) in our determination of flux. Since more 
concentrated halos are closer to the GC in our approach, 
the enhancement is greater close to the GC: indeed, at the 
anticenter it goes down to a factor 2. 

We have used the P(v(M)) for the ellipsoidal collapse in 
Eq. 1171 We have checked that using the corresponding prob- 
ability function for the spherical collapse does not change 
the result on $ cosmo . This is due to the fact that the main 
difference between the two models resides at small values of 
v. A small v gives low concentration parameter and its con- 




Figure 5. Subhalo contribution to the 7-ray flux for the two 
different models for the concentration parameters described in 
the text. MW smooth and clumpy contributions are shown sepa- 
rately, together with their sum. In the small box, zoomed at small 
angles from the Galactic Center only the sum is shown, and it is 
compared with the values obtained in lPieri et al.l (|2007| ). 



tribution to Eq.[T§]is then depressed with respect to that of 
a haloes with a higher v. 



4.3 Normalization to EGRET data 

In order to make predictions on detectability, we impose 
the best value of $ pp co mpatible w it h th e available 
experimental limits. As in IPieri et al.l ((2007), we first 
assume the optimistic model where m x = 40 GeV, 
CTaanV — 3 x 10 -26 cm 3 s _1 and the branching ra- 
tio is 100% in bb. We then integrate Eq |16l above 
3 GeV. This choice of parameters gives a value of 



$ pp = 2.6 x 10~ 9 cm 4 kpc" 1 GeV" 
We then compute the expected number of photons above 3 
GeV in 1 year for a solid angle of 10 -5 sr corresponding to 
the angular resolution of a GLAST-like satellite. The result 
for the B re /,o (dashed curve) and B re / :Z (dotted) models is 
shown in Fig. [6] 

We compare the obtained number of events with 
the EGRET data for t he diffuse Galactic co mponent 
parametrized according to iBergstrom et al l ((1998! ) 



dVtdE 



= N {l,b) 1(T 



7 



cm 2 s sr GeV ' 



(19) 



and with the diffuse extragalactic 7 emiss ion, as extrapo- 
lated from EGRET data at lower energies (|Sreekumar et all 
Il998t) : 



/ extra— 7 
diffuse 



1.38 x lO^iT 2 ' 1 7 - (20) 

dQ.dE cm 2 s sr GeV 

The normalization factor No in Eq. [T3] depends only on 
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Figure 6. Number of photons above 3 GeV, in 1 year in a solid 
angle of 10 -5 sr. The annihilation signal models B re ^o (dashed) 
and B re ^ z (dotted) are shown together with the EGRET diffuse 
expected Galactic and extragalactic background (solid) , as a func- 
tion of the angle of view ip from the Galactic Center. 



Figure 7. Sensitivity curves for a GLAST-like experiment, for 
the B re y o (solid) and the B re j> z (dotted) models described in the 
text. A zoom at small angles is provided in the superimposed 
frame. 



the interstellar matter distribution. The resulting number 
of photons above 3 GeV in 1 year for Afl = 10~ a sr, com- 
puted along 1=0 where its value is minimum, is shown in 
Fig. [6] (solid curve). 

We find an excess of annihilation signal photons toward the 
Galactic centre in both models. Yet, the angular resolution 
of EGRET corresponding to AQ = 10 -3 sr does not allow 
to reconstruct a spiky source as it is ours. We have checked 
that, if we compute the number of annihilation signal pho- 
tons toward ip = smeared in a cone of view of 1° , it is below 
the number of EGRET detected photons for the same an- 
gular resolution. 

Yet, the B re / iZ model exceeds the extragalactic diffuse mea- 
sured background too, which is dominant above tp = 40°. 
Since the extragalatic background is not due to any point 
source, we safely expect that it will scale with the solid angle. 
The number of annihilation signal photons produced in the 
B re /, z model should then be less or at most comparable with 
the number of measured background photons. We make the 
optimistic assumption that the two numbers are compara- 
ble at ip = 40° where the discrepance is larger, and we thus 
fix $If e/z = 2.0 x 10~ 9 cm 4 kpc- 1 GeV- 2 s- 1 sr- 1 for the 
B re / jZ model, correctly normalized to EGRET data, while 
we keep = 2.6 x KT 9 cm 4 kpc" 1 GeV -2 s" 1 sr" 1 for 

the B re /,o model. 



5 PROSPECTS FOR DETECTION 

In this section we study the sensitivity of a GLAST-like 
apparatus for 1 year of effective data taking. 
We define the experimental sensitivity a as the ratio of the 
number n 7 of annihilation signal photons and the fluctuation 



of background events 7ibk g : 



=- / Af(E, e t )[dcj>^ snal /dEdn]dEdQ. 

y/f E bkg KU E > Oi)ld<l> bkg /dEdCi\dEdn 

where Tg = 1 year is the effective observation time and <f>bkg 
is the background flux given by Ecis ll9l and 1201 computed 
along 1=0, that we assume to be composed by astrophysical 
photons only. The quantity eaq is the fraction of signal 
events within the optimal solid angle Af2 corresponding 
to the angular resolution of the instrument and it is 
optimistically set to 1. A eS is the effective detection area 
defined as the detection efficiency times the geometrical de- 
tection area. We use A eS = 10 4 cm , independent from the 
energy E and the incidence angle 9i. Finally we assume an 
angular resolution of 0.1° and an energy threshold of 3 GeV. 

The resulting sensitivity curves as a function of the an- 
gle of view ip are shown in FigQfor the H r ef,o (solid curve) 
and B re / iZ (dotted) annihilation signal models. In the small 
box a zoom at GC is shown. An almost 2 a around 10° is 
found for the B re / iZ model. The same model would be de- 
tected at about 30 a at the Galactic Centre. As far as the 
B re /,o model is concerned, it would show up with ~ 40cr 
effect toward the GC, that would rapidly fall down 1 a after 
0.5°. A 5 a detection at the Galactic Center would be possi- 
ble for both models with a value of <I> PP even 6 times lower. 
In case of a striking excess detection along the GC, a milder 
excess a larger angles could be a hint for the discrimination 
abo ut the models, tho ugh no discovery could be claimed. 

IPieri et all (|2007h studied the detectability of resolved 
haloes which would shine above the Galactic foreground, 
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finding in their best case scenario that only a tenth of large 
mass haloes would be detected, with a mass slope of -2 for 
the halo mass function. 

Repeating their analysis is beyond the goal of this pa- 
per. Yet we note that the effect of including the P(y(M)) 
factor in Eq.ll7lw ith respect to the concentration models in 
IPieri et al.l (|2007l ) leads to an enhancement of the Galactic 
foreground. We thus expect that including P(u(M)) will be 
compensated by the increased foreground and we don't ex- 
pect a dramatic change in the number of detectable haloes. 
As a further test, we have computed the sensitivity of 
a GLAST-like experiment for a B re / jZ halo once $pp has 
been normalized to the EGRET data. We chose the closer 
M = 1O~ 6 M halo, located at 9.5 x 1CT 2 pc from the sun. 
We chose v = 2.4 given from the probability of finding 1 
halo with such a value in a 1 pc 3 sphere around the sun. We 
conservatively considered only the astrophysical background 
in Eq. 1221 while the annihilation signal foreground should be 
considered too. Even in these very optimistic hypothesis, we 
found that the source would produce a 5 a effect only fur- 
therly multiplying by a factor of 4 the concentration param- 
eter. This could be achieved using the lognormal probability 
P(c(M)) but with a ridicolously small probability. 
We conclude that the effect of introducing the P(u(M)) 
can only be observed in a global enhancement of the dif- 
fuse Galactic annihilation foreground. 



6 CONCLUSIONS 

In this paper we have, for the first time, derived an an- 
alytical description of the mass function and distribution 
of rareness of density peaks in the subhalo population of 
our Galaxy, applying the extended Press & Schecter for- 
malism. To make the calculation possible, tidal interactions 
and close encounters between subhaloes have been neglected. 
Very small (micro solar mass) subhaloes are extremely con- 
centrated, therefore, at least for them, our approximation is 
a reasonable one. 

The obtained results are valid over the whole range of 
subhalo masses [10 -6 , 10 10 ]Mq and thus confirm and extend 
the results of the N-body simulations, whose resolution is 
still far too low in order to simulate coherently this mass 
range. 

Making use of the results of iDiemand etaH <|2005bl ) 
on the distribution of different a-peak material inside our 
Galaxy, we have been able to shape and model the total 
expected annihilation 7-ray foreground, statistically taking 
into account the merging history of each progenitor. 

We have used the best case particle physics scenario 
to derive predictions for the detectability of such a signal 
with a GLAST-like experiment. We have shown how both 
the merging history and the intrinsic properties of the halo 
formation can contribute to an enhancement of the expected 
flux, by arising the inner concentration of subhalos. Yet the 
real concentration of the single subhalo today remains an 
open question. We use two models which result in very dif- 
ferent inner densities inside the haloes. In the first model we 
assume that the inner shells of the subhaloes remain frozen 
at the moment they enter the parent halo and thus compute 
the concentration parameter at the merging epoch, as it is 
derived in our calculations. Alternatively we assume that the 



subhaloes continue to evolve with redsh ift, and thus compute 
the halo properties today. We use the lBullock et alj (|200ll ) 
model for the concentration par ameter at z = 0, extrapo- 
lated at low masses. We refer to IPieri et all l|2007l) for the 
effect of using different models. 

Our results on detectability show that a detection would 
be possible and impressive toward the GC for both models. 
This detection would be mainly due to the spike in the MW 
halo at the GC. Unfortunately, a reliable modeling of the 
astrophysical background coming from the GC and of the 
effect of the central Super Massive Black Hole on the inner 
DM density profile are still poorly known. 

A 2-<7 effect would show up as well, around ~ 10° from 
the GC, only for the B re / !Z model. Though no discovery could 
be claimed for, this could be a significant hint for the ex- 
istence of such a population of subhaloes, and it would be 
propulsive for successive studies with upcoming experimen- 
tal technologies. 

A final note on the metodology. In the present work 
we derived the final subhalo mass function starting from 
all progenitor haloes at any redshift. We did so in order 
to directly compare our analytica l results to the results 
obtained bv lDiemand et aD (|2005bl ) using JV-Body simula- 
tions. However, the subhalo population should indeed be 
derived starting from the population of "satellite haloes" 
directly accreted by the proto -halo (also cal led main pro- 
genitor) at all previous times (|Tormenlll997h . since only a 
fraction of progenitors at redshift 2 merge directly with the 
main halo progenitor. Unfortunately, the mass function of 
satellite haloes cannot be obtained analytically: it requires 
Monte Car lo simulations of the merging history tree of halo 
formation (ISomerville fc Kolattlll999l ; Ivan den Boschll2002l ; 
Ivan den Bosch et al.ll2005l ). We are currently working on this 
issue (Giocoli et al, in prep.), and it will be interesting to 
compare the results obtained using the two methods. 
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